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Abstract. - We investigate the dynamics of a modified Burridge-Knopoff model by intro- 
ducing a dissipative term to mimic the bursts of acoustic emission (AE) from rock samples. 
The model explains many features of the statistics of AE signals observed in experiments such 
as the crossover in the exponent value from relatively small amplitude AE signals to larger 
regime, and their dependence on the pulling speed. Significantly, we find that the cumulative 
energy dissipated identified with acoustic emission can be used to predict a major slip event. 
We also find a data collapse of the acoustic activity for several major slip events describable 
by a universal stretched exponential with corrections in terms of time-to-failure. 



Predicting failure of materials is of interest in science and engineering (electrical break- 
down, fracture of laboratory samples to engineering structures) particularly so in seismology 
due to the enormous damage earthquakes can cause. Whether it is at a laboratory or geolog- 
ical scale, this amounts to identifying useful precursors at a statistically significant level. One 
important non-destructive tool in fracture studies is the acoustic emission (AE) technique as 
it is sensitive to the microstructural changes occurring in the sample. Insight into earthquake 
dynamics has been obtained through fracture studies of ( usually precut samples to mimic slip 
on preexisting tectonic faults) rock samples [1]. Such studies have established that there is a 
considerable overlap between AE and seismology as both are concerned about the generation 
and propagation of elastic waves. Quite early, the statistics of the AE signals was shown to 
exhibit a power law [1-3] similar to the Gutenberg-Richters law for the magnitudes of earth- 
quakes [4] and Omori's law for aftershocks [2,5]. These prompted further investigations to 
look for precursor effects that can be used for earthquake predictability [5-8] . 

Apart from the power law distributions observed in the AE signals during fracture, acoustic 
activity of unusually large number of situations as varied as volcanic activity [9], micro- 
fracturing process [10,11], and collective dislocation motion [12], exhibits power laws statistics. 
Though the general mechanism attributed to AE is the release of stored strain energy, the 
details are system specific. Thus, the ubiquity of the power law statistics of AE signals 
suggests that the details of the underlying processes are irrelevant. Indeed, the conceptual 
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framework of self-organized criticality (SOC) was introduced to explain the universality of 
power law statistics in varied systems [13,14]. This brings up a related puzzle that lack of 
intrinsic length scales and time scales in power law systems has been taken to mean lack of 
predictability of individual avalanches [5, 15]. This generally held opinion has also triggered 
considerable debate regarding earthquake predictability, as seismically active fault systems 
are considered to be in a SOC state [15]. Thus, it would be desirable to design a model which 
exhibits power law statistics where predictability of individual events is possible. 

On the other hand, there has been reports of increased levels of AE signals before failure of 
rock samples [6, 16, 17] , in experiments on laboratory samples [18, 19] and precursory effects in 
individual earthquakes as well [5,7,20,21]. There has been attempts to predict failure within 
the framework of time-to-failure (TTF) models, both on the laboratory level [18] as well on 
geological scale [22]. (Some efforts have also been made on the predictability of avalanches 
in SOC state [23].) However, as the primary aim of TTF models is to mimic the rate of 
increase in precursory variable, they are independent of the nature of the variable. Thus, if 
one is interested explaining the origin of increasing level of activity in AE signals observed 
in experiments as the system approaches failure, it is necessary to first model AE in terms 
of displacement related variable as AE corresponds to high frequency elastic waves. Here, 
we model AE signals by introducing an additional dissipativc term into Burridge-Knopoff 
model [24-26] that helps us to identify precursory effect and hence to predict slip events. 

Our second objective is to explain some recent results on AE studies on rock samples. 
These studies report an interesting crossover in the exponent value from small amplitude 
regime to large [28] , a result that is similar to the well noted observation on the change in the 
power law exponent for small and large magnitudes earthquakes [29]. The exponent value is 
also found to be sensitive to the deformation rate [30]. To the best of our knowledge, there 
has been no explanation of these observations. This can partly be traced to the lack of efforts 
to model AE signals in terms of displacement related variables. 

Deformation and/or breaking of the asperities results in an accelerated motion of the local 
areas of slip where the stored potential energy is converted to kinetic energy. The general 
mechanism attributed to acoustic emission is the sudden release of the stored strain energy. In 
the case of modeling AE signals involving abrupt dislocation motion, the energy of AE signals, 
E ae (r) is taken to be proportional e 2 (r), where e is the local plastic strain rate [31]. However, 
in general there is spatial inhomogencity which implies that total energy E ae oc J (Ve) 2 d 3 r. 
For the BK model, displacement rate takes the role of e. On the other hand, E ae , has the same 
form as the Rayleigh dissipation functional [32] arising from the rapid movement of a localized 
region. Such a rapid movement prevents the system from attaining a quasi-static equilibrium 
which in turn generates dissipative forces that resist the rapid motion. On general grounds 
this dissipative forces has been shown to be proportional to / (\7v) 2 d 3 r, where v is the local 
velocity. Such a term is termed as solid viscosity in parallel with shear viscosity of fluids. 
(We ignore other types of dissipation such as radiation damping.) Indeed, we have found 
that the above term mimics all the features of the AE signals observed during martensitic 
transformation [33] ( both the power law statistics of the AE signals during thermal cycling 
and the highly correlated AE signals associated with growth and shrinkage of martensitc 
domains). 

The Burridge-Knopoff (BK) model for earthquakes [24] and its variants are popular models 
among the physics community. Despite its limitation (lack of appropriate continuum limit) 
[26,27], it forms a convenient platform to investigate the question of predictability of large 
avalanches as the dissipated energy bursts (in our modified model) themselves follow a power 
law distribution. The model consists of a chain of blocks of mass m coupled to each other by 
coil springs of strength k c and attached to a fixed surface by leaf springs of strength k p as 
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shown in Fig. ^ The blocks are in contact with a rough surface moving at constant speed v 
(mimicking the points of contact between two tectonic plates). A crucial input into the model 
is the velocity-dependent frictional force between the blocks and the surface (see Fig. 
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Fig. 1 - The Burridge-Knopoff spring block model and the two forms of friction laws used. 
Fig. 2 - Dissipative functional, R(t) as a function of t. 



The additional dissipative force is introduced through the Lagrange's equations of motion 
given by ^ ( sl(x) ) ~ Sufx) = ~ siifx) w ^h the Lagrangian L = T — P. The kinetic energy T is 
defined by T = f j(^-) 2 dx, and the potential energy by P = \ J[k c (^-) 2 + k p (u(x)) 2 }dx. 
Here u is the displacement of the blocks measured from the initial equilibrium position. The 
Rayleigh dissipative functional [32] is given by R = ^ /( 8 g^ ) 2 dx, where 7 is a dissipation 
coefficient. The total dissipation F is the sum of R(t) and frictional dissipation given by 
Ff r = j[f(v + ii(x))]dx, where the over dot refers to the time derivative and the frictional 
force is taken to be derivable from a potential like function. Then the equation of motion is 
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The discretized version, in the notation of Ref [25], reads 



U, 
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2Uj + Uj-i) - Uj - <p(2av + 2a(j j ) + j c (U j+ i - 2Uj + U^ x ), 



(1) 



(2) 



where Uj is the dimensionless displacement of the j block, v the dimensionless pulling 
velocity^ 2 = k c /k p the ratio of the slipping time to the loading time, and a is the rate of 
velocity- weakening in the scaled frictional force <fi (see Fig. The Coulomb frictional law 
is shown by the solid line. The dashed curve uses a resistive creep branch ending at v a 
(~ 1CF 7 here) beyond which the velocity weakening law operates. -f c is the scaled dissipation 
coefficient. The continuum limit of Eq. (2) exists for the creep branch (even in the absence 
of R(t)) which ensures a length scale below which all perturbations are damped [25]. Such 
a length scale is absent for the Coulomb case even as the continuum limit exists due to the 
additional dissipative term. These differences between the frictional laws may influence the 
nature of the statistics which we investigate in detail. 

This model without the last term has been extensively studied [24-26]. Starting from 
random initial conditions for all the blocks, slip events ranging from one-block event to those 
extending over the entire fault (occurring roughly once in a loading period T£ ~ 2/v ) are 
seen in the steady state for both types of frictional laws. 

Equation (2) has been solved using a fourth-order Runge-Kutta method with open bound- 
ary condition for both types of frictional laws shown in Fig. ^ Random initial conditions are 



4 



imposed. After discarding the initial transients, long data sets are recorded when the system 
has reached a stationary state. The parameters used here are I = 10, a = 2.5, TV = 100, 200 
for v = 0.01 and 0.001 and a range of values of j c . The modified BK model produces the 
same statistics of slip events as that without the last term in Eq. (2) as long as the value of 
7 C is small, typically 7 C < 0.5. The results presented here are for N = 100 and -f c = 0.02. 

Since the rate of energy dissipated [32] due to local accelerating blocks is given by dE ae /dt = 
—2R(t), we calculate R(t) which exhibits bursts similar to the AE signals observed in exper- 
iments. A plot of R(t) for the case with Coulomb frictional law is shown in Fig. [2] For the 
case with creep branch, R(t) is less noisy. 

We now consider the statistics of the energy bursts R(t). Denoting A to be the amplitude 
of R{t) (i.e., from a maximum to the next minimum), we find that the distribution of the 
magnitudes D(A) (accumulated over 2 x 10 5 time steps), shows a power law D(A) ~ A~ m . 
Instead of a single power law anticipated, we find that the distribution shows two regions, 
one for relatively smaller amplitudes and another for large values shown by the two distinct 
plots for the Coulomb case shown in Fig. ( The crossover in the power law exponent from 
smaller to larger amplitude region is seen as a scatter at the upper end of A ~ 10~ 5 .) Similar 
results are obtained when the frictional law has a creep branch. The value of m for the small 
amplitudes region ( typically < 10~ 4 ) is ~ 1.78±0.01, significantly smaller than that for large 
amplitudes which is ~ 2.09 ± 0.02, consistent with the recent experimental result on AE of 
rock samples [28] . The increase in the exponent value mimics a similar observation for large 
magnitude earthquakes ( > 7.0 on the Richter scale [29]). 

Recently Yabe et al. [30] have reported that the exponent value corresponding to relatively 
small amplitude regime increases with decreasing deformation rate while that for the large 
amplitude regime is found to be insensitive. To check this, we have calculated D(A) for 
v = 0.001 shown in Fig. |3 which shows that the exponent for small amplitude increases to 
1.91 ±0.02. However, we find that the exponent for larger amplitude regime is insensitive (not 
shown) to the changes in v. 

This result can be physically explained by analyzing the influence of the pulling velocity on 
slip events of varying sizes. We first note that R(t) depends on the difference in the velocities 
of neighboring blocks. The velocity of 'microscopic' events (small number of blocks) has been 
shown to proportional to v [25]. For single block events, as the neighboring blocks are at rest, 
the number of such events are fewer in proportion to the pulling speed, both of which are 
evident from Fig. |3| For the two block events, the contribution comes mostly from the edges 
as the difference in the velocities of the two blocks are of similar magnitude. In a similar way, 
it can be argued that for slip events of finite size, the extent of the contribution to R(t) is 
decided by the ruggedness of the velocity profile within the slipping region; the magnitude 
of R(t) is lower if the velocity is smoother. The ruggedness of the velocity profile, however, 
is itself decided by how much time the system gets to 'relax'. At lower values of v, there is 
sufficient time for the blocks to attain nearly the same velocities as the neighboring blocks 
compared to that at higher v values. Thus, the larger slip events contribute lesser to R(t) for 
smaller v values and hence the slope of logD(A) increases for smaller value of v. 

Now we consider the possibility of a precursor effect. A plot of R(t) is shown in the inset 
of Fig. 01 for the frictional law with a creep branch. A gradual increase in the activity of 
the energy dissipated can be seen which accelerates just prior to the occurrence of a 'major' 
slip event. The rapid increase in R(t) coincides with the abrupt increase in the mean kinetic 
energy (KE). Here, we have used KE as a measure of event size as it is a good indicator of 
the magnitude of the slip events. We find similar increase in R(t) for all 'major' slip events. 
(In our simulations, the KE of observable events ranges from 10~ 4 to 0.2. We refer to all such 
events as 'major' events.) This suggests that R(t) can be used as a precursor for the onset of a 
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Fig. 3 

Fig. 3 - Distribution of amplitudes of R(A), D(A) versus A. for small amplitudes (<^ for v — 0.01 • 
for v = 0.001 ) and large amplitudes (o for v — 0.01 shifted up for clarity). 

Fig. 4 - A semi- log plot of the cumulative energy, E ae versus t (o) along with the fit (solid line). 
Dashed line corresponds to the mean kinetic energy. The inset shows R(t) as a function of t. 



major slip event. As R(t) is noisy, a better quantity for the analysis is the cumulative energy 
dissipated E ae (t) (oc J Q R(t')dt'). E ae (t) grows in a stepped manner with their magnitudes 
increasing as we approach a major event. A plot of logE ae is shown in Fig. 01 along with a fit 
(continuous curve) having the functional form, 

logE ae {t) = -a x t~^ [l-a 2 \(t- t c )/t c \- a *]. (3) 

Here, t is the time measured from some initial point after a major event. The constants 
ai, 0,2, e*i and c*2 and t c are adjustable. The crucial parameter t c is the time of occurrence of 
a major slip event (often referred to as the 'failure' point). It is clear that the fit is striking. 
Given a reasonable stretch of the data, the initial increasing trend in logE ae is easily fitted 
to a stretched exponential, i.e., — ait~ ai . The additional term is introduced to account for 
the observed rapid increase in the activity as we approach the major event. In contrast, the 
mean kinetic energy increases abruptly on reaching the event ( dashed line). It is clear that 
the estimated t c agrees quite well with that of KE of the event. 

Now, we address the question of predictability of major slip events using E ae (t). This 
is equivalent to determining the correct t c . Given E ae (t) over a reasonable initial stretch of 
time, say till t\ (the first arrow in Fig. EJ, we find that the four constants ai,a,2,ai and a 2 
are already well determined (within a small error bar). These change very little with time. 

(Only t c changes.) A fit to Eq. (3) also gives tc at t\ ( which can only be considered as an 
estimate based on the data till t\). One such curve is shown by a dashed line with the arrow 
shown at t\. However, as time progresses, the data accumulated later usually deviates from 
the predicted curve if t c is inaccurate as is the case for the fits till t\ and t 2 for instance (Fig. 
EJ. If on the other hand, the deviation of the predicted curve from the accumulated data 
decreases with passage of time within the error bar (as is the case for the region just before 
ts), then, the value of t c is likely to be accurate. Indeed, the extrapolated continuous curve 
corresponding to data fit till t — t 3 (third arrow in Fig. 01 with the predicted is seen to 
follow the data very well. (Usually, the data deviates from the predicted curve with a sudden 
decrease in E~<t which is again an indication of a coherent slipping of several blocks before the 
onset of a fully delocalized event. But the general trend soon follows the extrapolated curve.) 
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Fig. 5 - A plot of logE a c versus t. Inset shows the enlarged section at time ti ( — ), %i ( — •) and 
(solid line). Data shown (•) is indistinguishable from the fit till t$. 

Fig. 6 - Collapsed data using aY 1 log(E ae (0)/ E ae (r)) vs. r for three different events along with the 
fit shown by (o). 

Then, £3 can be taken as the warning time for the onset of the major event. The actual t c 
read off from the kinetic energy plots is 111.0 where as the predicted t c is 111.6 giving the 
accuracy in the prediction of correct t c to be ~ 99.5%. This fit is obtained when ti as t is 12% 
away from the true t c . Similar numbers are obtained for several other events fitted. If the 
approach to all slip events is described by the scale invariant form, then one should expect to 
find a data collapse for different events. Indeed, in terms of a scaled time t = t/t Cl we find 
that the data corresponding to different events collapses into a single curve given by 

a^log[E ae {0)/E ae {T)} = [1 - a 2 |(r - 1)|- Q2 ] + a^logE ae (0). (4) 

A plot corresponding to three different events is shown in Fig. [fj] along with the fit. The 
results are similar when the Coulomb frictional law is used except that R(t) is more noisy and 
hence prone to slightly larger errors in predicted t c . 

In summary, the model captures several experimentally observed features of acoustic emis- 
sion such as the two different exponent values in the power law statistics for small and large 
amplitude regimes with the former being more sensitive to the pulling speed. The dependence 
of the exponent on the pulling speed has been traced to the form of R(t), namely, the gradient 
of the local velocity. More significantly, the analysis shows that it is possible to predict a 
major event fairly accurately. This result demonstrates that power law statistics in itself does 
not preclude predictability (of individual events). Further, we note that a SOC state demands 
that all observable quantities should follow a scale invariant form which is clearly respected by 
Eqs. (3,4) representing the approach to all events. The data collapse for different events into 
a scale free form clearly suggests that the dynamics of approach to major events is universal. 

A few comments are in order. We stress that this precursor effect is absent in the total 
kinetic energy or seismic moments. As there is no correspondence between seismic moments 
or the KE of the events with R(t) (which depends on the difference between the velocities of 
neighboring elements), our analysis cannot predict the magnitudes of the slip events. Indeed, 
in some cases we find that R(t) is larger for a smaller slip event, as is the case for the two 
peaks in R(t) shown in Fig. We point out here that the gradual increase in the energy 
dissipated as we approach a major event is different from that reported in the context of the 
BK model [34]. Our work also differs from the approach of Huang et al. [22] in the sense that in 
their analysis, the hierarchical structure and long range interaction are necessary ingredients 
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for the power law approach to failure with log periodic corrections that arises due to discrete 
scale invariance [19,22]. As far as we know, this is the first time that reproduces several 
unexplained experimental results on AE of rock samples [28,30]. We believe that it would be 
possible to address the problem of predictability in other power situations, which however, 
may demand capturing the essential physics of the relevant precursor effect. 
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work is supported by Department of Science and Technology, Grant No. Sp/s2K-26/98, New 
Delhi, India. 
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